library(foreign)
library(Zelig)
library(car)
library(mgcv)
library(VGAM)
library(ZeligChoice)
library(ordinal)
library(MASS)
library(stargazer)

setwd("C:/Users/Chan Peter Kim/Desktop/replication paper")

kroenig <- read.dta("Nuclear Crisis Outcomes updated.dta")
kroenig <- kroenig[1:52,]


############ TABLE 1 ###############

clx <- 
  function(fm, dfcw, cluster){
    library(sandwich);library(lmtest)
    M <- length(unique(cluster))   
    N <- length(cluster)           
    K <- fm$rank                        
    dfc <- (M/(M-1))*((N-1)/(N-K))
    uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
    vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)*dfcw
    coeftest(fm, vcovCL)
    return(vcovCL)
  }


model1.1.formula <- as.formula(victory ~ nuc_rat + proximity + stakes + polity21
                               + capshare + strikea + tpop_1 + viol + cris_ave)
model1.1 <- glm(model1.1.formula, data = kroenig, family=binomial(link="probit"))
summary(model1.1)$coefficients
model1.1.se <- clx(model1.1, 1, kroenig$crdynum)
model1.1.se <- sqrt(diag(model1.1.se))


model1.2.formula <- as.formula(factor(victory_ordinal) ~ nuc_rat + proximity + stakes + polity21
           + capshare + strikea + tpop_1 + viol + cris_ave)
model1.2 <- polr(model1.2.formula, data = kroenig, Hess = T, method = "probit")
model1.2.coef <- model1.2$coefficients
sqrt(diag(model1.2$Hess))
ctable <- coef(summary(model1.2))
model1.2.pvals <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
model1.2.pvals


################## TABLE 2  ################

model2.1.formula <- as.formula(factor(victory_ordinal) ~ nuc_rat + proximity + stakes + polity21
                               + milex_rat + strikea + tpop_1_new + viol
                               + cris_ave + electionyear)
model2.1 <- polr(model2.1.formula, data = kroenig, Hess = T, method = "probit")
model2.1.coef <- model2.1$coefficients
sqrt(diag(model2.1$Hess))
ctable <- coef(summary(model2.1))
model2.1.pvals <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
model2.1.pvals

model2.2.formula <- as.formula(factor(victory_ordinal) ~ superiority + proximity + 
                                 stakes + polity21 + milex_rat + strikea + tpop_1_new + viol
                               + cris_ave + electionyear)
model2.2 <- polr(model2.2.formula, data = kroenig, Hess = T, method = "probit")
model2.2.coef <- model2.2$coefficients
sqrt(diag(model2.2$Hess))
ctable <- coef(summary(model2.2))
model2.2.pvals <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
model2.2.pvals

stargazer(model2.1, model2.2, ord.intercepts = T, font.size = "small")



####### DURATION MODELS ###########


library(xtable)
library(MASS)


ll.dur <- function(par, y, c, X){
  beta <- par[1:ncol(X)]
  lambda <- exp(X%*%beta)
  alpha <-  exp(par[ncol(X)+1])
  sum((1-c)*(log(alpha) - alpha*log(lambda) + (alpha-1)*log(y) - (y/lambda)^alpha)
      - c*(y/lambda)^alpha)
}


## ASYMMETRY ##

X <- cbind(1, as.matrix(kroenig[,c("asymmetry2", "proximity", "stakes", "polity21", "milex_rat", "strikea", "viol", "cris_ave", "nato", "demvsnondem")]))
c <- as.matrix(rep(0, nrow(kroenig)))
y <- as.matrix(kroenig[,"durdays_new"])

my.opt <- optim(par = rep(1, ncol(X) + 1), fn = ll.dur, y = y, c=c,
                X = X, control = list(fnscale = -1), method = "BFGS", hessian = TRUE)

my.opt.coefs <- my.opt$par
my.opt.ses <- sqrt(diag(solve(-my.opt$hessian)))

res <- cbind(my.opt.coefs, my.opt.ses)
colnames(res) <- c("Coefficient","SE")
rownames(res) <- c("Intercept", "asymmetry2", "proximity", "stakes", "polity21", "milex_rat", "strikea", "viol", "cris_ave", "nato", "demvsnondem", "Alpha")
xtable(res, digits=4)

# Simulate Parameters
set.seed(01238)
simpar <- mvrnorm(10000, my.opt$par, -solve(my.opt$hessian))

# Set Covariates at mean
setx <- cbind(apply(X,2,median))

# Simulate Fundamental Uncertainty
sims <- NA
for(i in 1:nrow(simpar)){sims[i] <- rweibull(1, shape=exp(simpar[i,12]), scale=exp(simpar[i,1:11]%*%setx))}

# Set Covariates 
setx.asymmetry <- c(1, 1, apply(X[,3:11], 2, median))
setx.symmetry <- c(1, 0, apply(X[,3:11], 2, median))

# Simulate Fundamental Uncertainty and Calculate expected values
set.seed(02138)
evs <- matrix(NA, nrow=nrow(simpar), ncol=2)
for(i in 1:nrow(simpar)){
  simfund <- rweibull(1000, shape=exp(simpar[i,12]), scale=exp(simpar[i,1:11]%*%setx.asymmetry))
  evs[i,1] <- mean(simfund)
  simfund2 <- rweibull(1000, shape=exp(simpar[i,12]), scale=exp(simpar[i,1:11]%*%setx.symmetry))
  evs[i,2] <- mean(simfund2)
}

mean(evs[,1]) - mean(evs[,2])

pdf("duration_asymmetry.pdf")
plot(density(evs[,2]), main="Nuclear Asymmetry", xlab="Expected Duration of Nuclear Crisis", col="orangered1", lwd=3, xlim=c(0,200), ylim = c(0, .02), cex.main = 2, cex.axis = 1.2, cex.lab = 1.4, cex.sub = 1.5)
lines(density(evs[,1]), col="royalblue3", lwd=3)
legend(x="topright", legend=c("Asymmetry", "Symmetry"), col=c("royalblue3", "orangered1"), lty=1, lwd=3, cex=1.3)
abline(v = mean(evs[,2]), col = "orangered1", lty = 2, lwd = 2)
abline(v = mean(evs[,1]), col = "royalblue2", lty = 2, lwd = 2)
dev.off()







### SECOND STRIKE ##

X <- cbind(1, as.matrix(kroenig[,c("strikebalance", "proximity", "stakes", "polity21", "milex_rat", "viol", "cris_ave", "nato", "demvsnondem")]))
c <- as.matrix(rep(0, nrow(kroenig)))
y <- as.matrix(kroenig[,"durdays_new"])

my.opt <- optim(par = rep(1, ncol(X) + 1), fn = ll.dur, y = y, c=c,
                X = X, control = list(fnscale = -1), method = "BFGS", hessian = TRUE)

my.opt.coefs <- my.opt$par
my.opt.ses <- sqrt(diag(solve(-my.opt$hessian)))

res <- cbind(my.opt.coefs, my.opt.ses)
colnames(res) <- c("Coefficient","SE")
rownames(res) <- c("Intercept", "strikebalance", "proximity", "stakes", "polity21", "milex_rat", "viol", "cris_ave", "nato", "demvsnondem", "Alpha")
xtable(res, digits=4)


# Simulate Parameters
set.seed(01238)
simpar <- mvrnorm(10000, my.opt$par, -solve(my.opt$hessian))

# Set Covariates at mean
setx <- cbind(apply(X,2,median))

# Simulate Fundamental Uncertainty
sims <- NA
for(i in 1:nrow(simpar)){sims[i] <- rweibull(1, shape=exp(simpar[i,11]), scale=exp(simpar[i,1:10]%*%setx))}

# Set Covariates 
setx.secondstrike <- c(1, 2, apply(X[,3:10], 2, median))
setx.mutualsecondstrike <- c(1, 1, apply(X[,3:10], 2, median))
setx.nosecondstrike <- c(1, 0, apply(X[,3:10], 2, median))

# Simulate Fundamental Uncertainty and Calculate expected values
set.seed(02138)
evs <- matrix(NA, nrow=nrow(simpar), ncol=3)
for(i in 1:nrow(simpar)){
  simfund <- rweibull(1000, shape=exp(simpar[i,11]), scale=exp(simpar[i,1:10]%*%setx.secondstrike))
  evs[i,1] <- mean(simfund)
  simfund2 <- rweibull(1000, shape=exp(simpar[i,11]), scale=exp(simpar[i,1:10]%*%setx.nosecondstrike))
  evs[i,2] <- mean(simfund2)
  simfund3 <- rweibull(1000, shape=exp(simpar[i,11]), scale=exp(simpar[i,1:10]%*%setx.mutualsecondstrike))
  evs[i,3] <- mean(simfund3)
}

mean(evs[,1]) - mean(evs[,2])

pdf("duration_secondstrike.pdf")
plot(density(evs[,2]), main="Second Strike", xlab="Expected Duration of Nuclear Crisis", col="orangered1", lwd=3, xlim=c(25,200), ylim = c(0, .02), cex.main = 2, cex.axis = 1.2, cex.lab = 1.4, cex.sub = 1.5)
lines(density(evs[,1]), col="royalblue3", lwd=3)
lines(density(evs[,3]), col="darkgoldenrod1", lwd=3)
legend(x="topright", legend=c("One State Second Strike", "Both States Second Strike", "Neither State Second Strike"), col=c("royalblue3", "darkgoldenrod1", "orangered1"), lty=1, lwd=3, cex=1.1)
abline(v = mean(evs[,2]), col = "orangered1", lty = 2, lwd = 2)
abline(v = mean(evs[,1]), col = "royalblue2", lty = 2, lwd = 2)
abline(v = mean(evs[,3]), col = "darkgoldenrod1", lty = 2, lwd = 2)
dev.off()